First-principles Simulation of Electrochemical Systems at Fixed Applied Voltage: 
Vibrational Stark Effect for CO on Platinum Electrodes 
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Chemisorbed molecules at a fuel-cell electrode are a very sensitive probe of the surrounding 
electrochemical environment, and one that can be accurately monitored with different spectroscopic 
techniques. We develop a comprehensive electrochemical model to study molecular chemisorption 
at fixed applied voltage, and calculate from first principles the voltage dependence of vibrational 
frequencies (the vibrational Stark effect) for CO on platinum electrodes, finding excellent agreement 
with electrochemical spectroscopic experiments and resolving previous controversies. 
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Rising sustainability concerns have revived strong in- 
terest in electrochemical electricity generation [1] whose 
basic principle is to catalytically convert the energy 
stored in chemical bonds into usable electrical power. For 
any given electrochemical system (e.g., a fuel cell or a 
battery), the power generated is the product of two dis- 
tinct contributions: (1) the electrode voltage difference, 
which is the thermodynamic variable that quantifies the 
energy per electron made available through the breaking 
and rearranging of chemical bonds (Nernst's law), and (2) 
the current density, which is the kinetic observable that 
measures the rate at which these chemical processes take 
place (Arrhenius' law). It should be emphasized, how- 
ever, that these two factors are not completely indepen- 
dent, as one observes experimentally a systematic drop 
in voltage at high electrical current. This phenomenon, 
commonly known as activation voltage loss, represents 
one of the main limitations to the performance of elec- 
trochemical technologies Q. 

Although the origins of the voltage dependence of the 
electrical current have long been conceptually under- 
stood, it is only recently that computational laboratories 
have applied first-principles calculations to study this ef- 
fect with the difficult task of describing catalytic reac- 
tions at an electrode surface as a function of the applied 
voltage. The electrochemical free-energy correction in- 
troduced by N0rskov et al. [!, Q represents a key success- 
ful step in this direction. In this approach, the influence 
of the electrode voltage £ is included by adding a correc- 
tion —e£ to the energy of all reaction intermediates that 
involve an electron transferred to the metal. This zeroth- 
order correction has been shown to accurately predict the 
activation voltage of fundamental electrocatalytic pro- 
cesses, such as the oxygen reduction reaction at fuel-cell 
cathodes However, this approach does not capture 
the self-consistent modifications of the electronic struc- 
ture that arise from the applied potential. In particular, 
it does not consider the variation of the electrode charge 
as a function of the potential and the interaction of the 



induced surface electric field with chemisorbed molecules. 

Several other authors have proposed to account for 
these important electronic effects using more represen- 
tative electrochemical models [5, 6|, 0, U 0] • These calcu- 
lations differ in key quantitative and qualitative details 
from the approach we present, and are generally carried 
out at constant electrode charge q, the voltage £ being 
determined a posteriori using various procedures to relate 
its value to the computed Fermi level ep. This is in clear 
contrast to most in-situ experiments, in which the state 
of the system is directly controlled via the electrode volt- 
age and all relevant electrochemical properties are given 
in terms of this central intensive variable. Although the 
voltage dependencies of electrochemical properties can be 
recovered via an inverse Legendre transform, this indirect 
method entails repeated constant-charge calculations to 
invert the charge-to- voltage relation. 

In this study, we introduce a practical computational 
model that allows to work directly at fixed electrode volt- 
age while fully describing self-consistent changes in the 
electronic structure and taking into account realistic elec- 
trochemical conditions. A validation of the method is 
provided by the prediction of the voltage dependence of 
vibrational frequencies — the vibrational Stark effect — for 
CO on platinum electrodes, which offers a very sensitive 
spectroscopic probe of the electrochemical environment 
and surface electric field. Since the vibrational proper- 
ties of chemisorbed molecules can be accurately described 
from first principles and precisely measured using 
various infrared techniques ll|, [l2[ , the calculation of the 



Stark effect represents a stringent test for assessing the 
predictive ability of an electrochemical model. To date, 
very large discrepancies between first-principles predic- 
tions and Stark measurements have been reported, and 
their elucidation has been a long-standing question in 
surface science and electrochemistry [l3j |. 

To put matters into perspective, a typical electrode- 
electrolyte interface is depicted in Fig. []Ja). The system 
consists of an adsorbate-covered metal surface in contact 
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FIG. 1: (a) Adsorbate-covered catalytic electrode-electrolyte 
interface and (b) implicit atom-continuum model of the dou- 
ble layer interface at fixed electrode voltage. 



with an electrolyte solution. For this system, the elec- 
trode voltage £ corresponds to the energy involved in 
displacing an electron from the metal electrode to the 
bulk of the ionic solvent |l4|]. Therefore, including volt- 
age conditions requires to accurately describe the behav- 
ior of the screened electrostatic potential in the solvent 
region. It is important to note, however, that electro- 
static screening in the electrolyte occurs on considerably 
large length scales — typically, 10-10 3 A for ionic con- 
centrations in the range 10 _4 -10 _1 mol/liter (M) [15[ — , 
which renders the explicit first-principles representation 
of the electrolyte prohibitively expensive and practi- 
cally inaccessible with current computational resources. 
Furthermore, quantitative errors in local and semilo- 
cal density-functional theory (DFT) descriptions of wa- 
ter [16[ have be reported. These errors result in over- 
structured representations of aqueous media and over- 
estimated freezing temperatures [17j — which translates 
into overevaluated dielectric responses for explicit solva- 
tion models. 

To overcome these important limitations, we introduce 
an atom-continuum double layer model of the electrified 
interface [Fig. [U(b)] [18[. This model consists of immers- 
ing the metal electrode (the explicit metal layer) in a 
semi-infinite electrolytic continuum (the implicit diffuse 
layer) [14|. Note that the thickness Ah of the interphase 
region (the Helmholtz interface) is experimentally found 
to be equal to 3-5 A (that is, approximately the thick- 
ness of a water bilayer) regardless of the nature of the 



electrode compound [15[. Within this implicit approach, 
the electrode voltage £ can be simply expressed as 
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where Vo = V is the asymptotic value of the poten- 
tial in the bulk of the electrolyte. With this physical 
picture in mind, it clearly appears that increasing the 
electrode voltage causes a depletion of surface electronic 
states compensated by an accumulation of negative coun- 
terions in the electrolyte [Fig. Gib)]. The induced po- 
larization enhances the double layer electric field, which 
interacts more strongly with the adsorbates and shifts 
their vibrational frequencies. 

In this model, a smoothly varying dielectric permit- 
tivity e accounts for the water environment, and diffuse 
charge densities and c~j~ represent the thermal dis- 
tributions of the counterions of bulk concentration Co, 
absolute charge z&, and size a& (for the sake of simplic- 
ity, we restrict ourselves to the case of a z&\z& symmetric 
ionic solution). The dielectric permittivity is calculated 
using the parametrization of Gygi and Fattebert that in- 
volves the static dielectric constant of water eo = 78 pjl • 
This dielectric model properly captures the gradual tran- 
sition of the permittivity across the solvation shell [2(j 
and yields accurate solvation energies for a broad range 
of molecular species 21]. The ionic concentrations c~j~ 



and c~,~ follow the modified Boltzmann statistics intro- 
duced by Borukhov, Andelman, and Orland [22[, which 
includes finite-size steric interactions between counteri- 
ons. Note that the contribution from explicit water over- 
layers can always be included. However, due the weak 
chemical interactions between CO and the solvent (23[, 
we restrict here the explicit DFT treatment to the metal 
and to the chemisorbed molecules. 

The ground state of the electrochemical system is that 
which minimizes the free energy functional 



G = E' + AE C 
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The first contribution E' corresponds to the DFT energy 
of the system within the supercell approximation (that 
is, for a periodically repeated metal slab in vacuum), as 
computed by standard plane- wave codes [24|. The cor- 
rective energy AE corr equals the difference between the 
electrostatic energy of the isolated slab and that of the 
periodic slab in vacuum, 



AE C 



p(r)?; corr (r)dr, 
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where p is the explicit charge density, and v corr = v — v' 
is a corrective potential defined as the difference between 
the Coulomb potential of the solvated system v that 
satisfies a nonlinear modified Poisson-Boltzmann (MPB) 
equation 

V • e(r)Vv(r) = -4n\p(r) - z d c+(r) + z d c^(r)} (4) 
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(in a.u.) with boundary conditions v = V at infinity, 
and the periodic potential v' calculated in the reciprocal 
space representation using fast Fourier transform tech- 
niques. Heuristically, v covv can be identified as the elec- 
trolyte reaction field [20 ]. which self- consistently accounts 
for the influence of the applied electrode voltage. The re- 
maining contribution AE lon is that from the counterions 
in solution, including steric repulsion [22| : 



AE 1( 
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-T J dr [s ion (c+, c", c° d ) - s ion (c , c , a" 3 )] , (5) 



where \i = — /cBTln(a^~ —2) is the ionic potential, s lon 
is the local ionic entropy, and is the saturated (max- 
imum packing) ionic concentration that smoothly goes 
from to a^ 3 at a distance Ah from the metal surface. 
The entropy density s lon can be expressed as 



s ioa (c+, c d -,c5) = -fc B 
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and the maximal concentration cj(r) is parametrized as 

<%<*)= \ n /§ (i r - R ^- A H) 5 ( 7 ) 

where the index / runs exclusively over the metal layer 
atoms (at position R/), and the counterion exclusion re- 
gion is defined by a smooth step function 6 (smeared 
over a few grid points for numerical convergence). Con- 
sequently, the equilibrium ionic concentrations read 
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{ Co - 1 a- 3 + 2[cosh(^)-l]}" 1 .(8) 



This ionic model directly involves the Helmholtz thick- 
ness Ah through the prefactor c£ and provides a simple 
representation of the diffuse distributions in direct con- 
nection to the Stern picture 

The implementation of this electrochemical model 
raises three main difficulties. First, solving the MPB 
electrostatic problem in a finite simulation cell with, e.g., 
periodic or homogeneous boundary conditions results in 
significant errors in the electrode voltage. To correctly 
extrapolate the slowly vanishing electrostatic potential, 
we choose to impose fictitious electrochemical boundary 
conditions obtained from the long-range integration of 
the MPB equation in the planar- average approximation: 
Vv • n z = _ / 327rc fc B T sinh ^v_ ( where Uz denotes the 



MPB equation is expensive due to the fine grids required 
in discretizing the charge density. To reduce this compu- 
tational burden, we exploit the fact that the corrective 
potential v corr varies smoothly over space, which allows 
for its direct and inexpensive computation on coarse- 
grained meshes in the spirit of the density-countercharge 
method [25[. Last, constant-potential simulations require 
fixing the Fermi energy while readjusting the electron 
number during the electronic-structure optimization. In 
the course of such calculations, large charge oscillations 
occur, which results in systematic energy divergence. To 
eliminate these instabilities without resorting to artificial 
charge compensations, we employ a generalization of the 
ensemble density-functional theory scheme [26] and opti- 
mal damping algorithm [27|, which ensures that the free 
energy converges monotonically. 
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external surface vector). Additionally, the solution of the 



FIG. 2: Vibrational frequency v as a function of the NHE- 
referenced electrode voltage 6 — 8° for 1/4 ML and 3/4 ML 
of CO on Pt(lll). The absolute electrode potential 8 is also 
indicated. The experimental Stark tuning rates dv/d8 are 
given in parentheses. 

We thus proceed to calculate the vibrational properties 
of CO-covered platinum interfaces under electrochemical 
conditions [3l|. In carrying out these calculations, the 
solvation parameters are set to the values determined 
and used in Ref. [2l|. The electrode voltage range, coun- 
terion concentration and ionic temperature are selected 
to be £ = 5.0-5.2 V, c = 0.1 M and T = 300 K, re- 
spectively, corresponding to experimental conditions. We 
reference the absolute potential of the half cell to that of 
the normal hydrogen electrode (NHE) by matching the 
electrode voltage calculated at the point of zero charge 
£pzc = 5.07 V to the referenced experimental potential 
£pzc P) - 8° = 0.33 V [28]. This corresponds to shifting 
the absolute electrode potential by 8° = 4.74 V in our 
calculations. The thickness of the double layer Ah equals 
4 A, that is, in the middle of the experimental range 3-5 
A. The size and charge of the counterions are chosen to be 
aa = 2 A and = 1 for a typical monovalent electrolyte. 
We compute stretching frequencies using a frozen-phonon 
method — the calculated frozen-phonon frequencies agree 
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to within 1-2 cm -1 to the results of the full computation 
and diagonalization of the dynamical matrix [10| . 

The dependency of the C-0 intramolecular frequencies 
as a function of the electrode voltage at surface concen- 
trations of 1/4 and 3/4 ML for the experimentally ob- 
served atop configuration is depicted in Fig. [2l First, 
we note that the predicted vibrational frequencies follow 
a common increasing and almost linear trend as a func- 
tion of the electrode potential, in qualitative agreement 
with experiment under preoxidation conditions (i.e., be- 
low ~0.5 V vs. NHE). In addition, the vibrational Stark 
effect is seen to strongly depend on the CO monolayer 
concentration. Indeed, at a coverage of 1/4 ML, the vi- 
brational Stark slope is calculated to be dv/dE = 43.2 
cm^/V, whereas at 3/4 ML, du/dS equals 28.4 cm _1 /V 
[321 ]. Despite this marked dependency as a function of 
monolayer coverage, the Stark rates are found to be in 
remarkable accordance with their experimental counter- 
parts du/dS = 40 cm" 1 /Vat 1/4 ML @ and 28 cm" 1 /V 
at 3/4 ML [HEH. Thus, at variance with previously re- 
ported large discrepancies, the calculated Stark tuning 
rates are found here to deviate by less than 8% from ex- 
perimental spectroscopic measurements in both the low- 
and high-coverage regimes, providing an important illus- 
tration of the predictive performance of the present elec- 
trochemical model. 

In summary, we have developed a practical and com- 
prehensive electrochemical model to study quantum- 
mechanical systems as a direct function of the fixed ap- 
plied voltage. We used this model to calculate voltage- 
induced Stark shifts, finding excellent agreement with 
spectroscopic experiments at low and high surface cov- 
erage, and resolving long-standing inconsistencies in the 
interpretation of electrochemical spectroscopic measure- 
ments. These results open promising perspectives for 
the first-principles description of electrochemical systems 
and electrocatalytic reactions under realistic voltage con- 
ditions. 

The calculations in this work have been performed us- 
ing the Quantum-Espresso package [30]. The authors 
acknowledge support from the MURI grant DA AD 19- 
03-1-0169 and Sire grant ANR 06-CIS6-014. We thank 
Andrzej Wieckowski for his help in providing and in in- 
terpreting spectroscopic data. Helpful suggestions and 
comments from Gerbrand Ceder, Stefano de Gironcoli, 
Damian Scherlis, Nicephore Bonnet, Jean-Sebastien Fil- 
hol, and Eduardo Lamas are gratefully acknowledged. 



* Electronic address: daboi@cermics. enpc.fr| 
[1] J. E. Tester et a/., Sustainable Energy: Choosing Among 
Options (MIT Press, 2005). 



[2] J. Larminie and A. Dicks, Fuel Cell Systems Explained 

(Wiley- VCH, 2003), 2nd ed. 
[3] J. K. N0rskov et al., J. Phys. Chem. B 108, 17886 (2004). 
[4] G. S. Karlberg et a/., Phys. Rev. Lett. 99, 12610 (2007). 
[5] A. Y. Lozovoi, A. Alavi, J. Kohanoff, and R. M. Lynden- 

Bell, J. Chem. Phys. 115, 1661 (2001). 
[6] A. Y. Lozovoi and A. Alavi, Phys. Rev. B 68, 245416 

(2003). 

[7] M. Otani and O. Sugino, Phys. Rev. B 73, 115407 (2006). 
[8] C. D. Taylor, S. A. Wasileski, J.-S. Filhol, and M. Neu- 

rock, Phys. Rev. B 73, 165402 (2006). 
[9] R. Jinnouchi and A. B. Anderson, Phys. Rev. B 77, 

245417 (2008). 

[10] I. Dabo, A. Wieckowski, and N. Marzari, J. Am. Chem. 

Soc. 129, 11045 (2007). 
[11] F. Maillard, G. Q. Lu, A. Wieckowski, and U. Stimming, 

J. Phys. Chem. B 109, 16230 (2005). 
[12] G. Q. Lu, A. Lagutchev, D. D. Dlott, and A. Wieckowski, 

Surf. Sci. 585, 3 (2005). 
[13] A. Y. Lozovoi and A. Alavi, J. Electroanal. Chem. 607, 

140 (2007). 

[14] J. O. Bockris and S. U. M. Khan, Surface Electrochem- 
istry (Plenum Press, 1993). 

[15] N. Sato, Electrochemistry at Metal and Semiconductor 
Electrodes (Elsevier, 1998). 

[16] J. G. Grossman et a/., J. Chem. Phys. 120, 300 (2004). 

[17] P. H.-L. Sit and N. Marzari, J. Chem. Phys. 122, 204510 
(2005). 

[18] I. Dabo, Ph.D. thesis, Massachusetts Institute of Tech- 
nology (2007). 

[19] J.-L. Fattebert and F. Gygi, J. Comput. Chem. 23, 662 
(2002). 

[20] A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. 

Chem. Theory Comput. 4, 877 (2008). 
[21] D. A. Scherlis et a/., J. Chem. Phys. 124, 74103 (2006). 
[22] I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. 

Lett. 79, 435 (1997). 
[23] A. Roudgar and A. Gross, Chem. Phys. Lett. 409, 157 

(2008). 

[24] M. C. Payne et a/., Rev. Mod. Phys. 64, 1045 (1992). 
[25] I. Dabo, B. Kozinsky, N. E. Singh-Miller, and N. Marzari, 

Phys. Rev. B 77, 115139 (2008). 
[26] N. Marzari, D. Vanderbilt, and M. C. Payne, Phys. Rev. 

Lett. 79, 1337 (1997). 
[27] E. Cances and C. Le Bris, Int. J. Quant. Chem. 79, 82 

(2000). 

[28] R. Gomez, V. Climent, J. M. Felium, and M. J. Weaver, 

J. Phys. Chem. B 104, 597 (2000). 
[29] M. J. Weaver, S. Zou, and C. Tang, J. Chem. Phys. Ill, 

368 (1999). 

[30] S. Baroni et a/., http://www.quantum-espresso.org/. 

[31] We employ the Perdew-Burke-Ernzerhof approximation 
and ultrasoft pseudopotentials. The Brillouin zone is 
sampled with a shifted 4 x 4 x 1 . Plane- wave energy cut- 
offs of 25 and 200 Ry. The system is represented by fully 
relaxed three-layer-thick Pt(lll) slabs at the calculated 
bulk lattice parameter of 3.99 A in a y/3 x 2 supercell. 

[32] By performing a sensitivity analysis, we find the Stark 
slopes to be altered by less than 3 cm -1 /V when varying 
the main experimental parameter Ah by ±1 A. 



